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Abstract 

In this paper, we presented an efficient algorithm to implement the regularization reconstruction of 
SPECT. Image reconstruction with priori assumptions is usually modeled as a constrained optimization 
problem. However, there is no efficient approaches to solve it due to the ill-posedness and large scale of the 
problem. In this paper, we use the superiorized approach of the expectation maximization (EM) iteration 
to implement the regularization reconstruction of SPECT. We first investigated the convergent conditions of 
the EM iteration in the presence of perturbations. Then we illustrated the superiorized EM algorithm based 
on the convergent conditions, and proposed a modified version of it. Furthermore, we gave two concrete 
methods to generate perturbations for two special objective functions. Numerical simulations for SPECT 
reconstruction were conducted to validate the performance of the proposed algorithms. The simulations show 
that the superiorized EM algorithms are more stable and robust for noise pollution and selection of initial 
point than the classic EM algorithm, and the reconstructed images by the proposed algorithms outperform 
the classic EM algorithm in terms mean square error and visual quality. 

Keywords: EM algorithm, superiorization, SPECT. 



1 Introduction 

Single-photon emission computed tomography (SPECT), which can visualize the physiological information of 
various organs with the help of radiopharmaceuticals [I] |21 [3] ■ The gamma-rays emitted by the injected radioac- 
tive material are recorded by a gamma camera rotated around the patient. The goal of SPECT is to reconstruct 
the radionuclide distribution from the measurements numerically. The variety of existing SPECT algorithms 
can be split into a family of analytical methods [H [5l |6j B |8] and a wide class of iterative techniques O [TOl [11] . 

From the analytic point of view, the SPECT reconstruction problem O HI [7] is to inverse the attenuated 
Radon transform (aRt) of /(distribution of radiopharmaceutical) 

JR 

where /i is a known function, referred to as the attenuation map of gamma-rays, 6 = (cos (f, sin ip) and 6-^ = 
(— sin(/5,cos<p). In practice, / is a function within compact support ft. Therefore, the integrand in (jlj is zero 
outside a bounded interval; this integral is written over (— oo, -l-oo) for convenience. 

For the iterative method, the SPECT reconstruction problem is modeled as the inversion of the following 
linear system [5], 

Ax = b, (2) 

where the elements of the observed data b = {bi,b2 ■ ■ ■ , bm) G M*^, the unknown image x = (a;i, X2 - ■ ■ , xn) & 
and the system matrix A = (aij) G ^mxn ^^^.^ nonnegative. The aim is to reconstruct the unknown x as 
an image from the projection data b via stable algorithms. A solution is not feasible with conventional methods 
directly because of the noisy projection data b, ill-posedness and large scale of the problem. 

For iterative methods, the system matrix A not only can model the attenuation of photon, but also can 
fuse some realistic factors, such as photon scattering and camera blurring. Therefore, we focus on the iterative 
methods in this paper. The expectation maximization (EM) algorithm [Sj [H] and the algebraic reconstruction 
technique (ART)[T31[I1] are two widely used technologies in imaging science, due to their simplicity, efficiency 



1 



and performance. 

In practice, the EM algorithm is more appropriate for emission tomography including SPECT. Firstly, the 
EM algorithm maintains the nonncgative constraint in the iteration procedure. Secondly, the EM algorithm 
is relatively robust against data inconsistencies introduced by Poisson noise, because it seeks to minimize the 
Kullback-Liebler(K-L) distance between the measured data b and the projection of the estimated image Ax, 
which is equivalent to maximize the likelihood of Poisson distribution. 

For SPECT reconstruction, it is one of the main issues to estimate the radionuclide distribution from low- 
counts projection data. This issue occurs quite frequently because of practical constraints, such as imaging 
hardware and scanning geometry. Furthermore, the decrease of the counts can reduce acquisition time, radia- 
tion dose and imaging cost efficiently. However, this would results the strong deterioration of the observed data 
and the under-determinacy (to ^ n) of the linear system ([2]). In these situations, the reconstructed images 
by the EM algorithm are usually dominated by various distortions [TJ [TSJ [121 E] ) because the EM algorithm 
accepts any solution which minimizes the K-L distance. 

The qualities of the reconstructions can be improved by regularization reconstruction methods. The reg- 
ularization reconstruction for SPECT is usually modeled as a large scale constrained optimization problem. 
However, to our knowledge, there is no optimal algorithms to solve this problem efficiently. In this paper, we 
resort to an emerging approach called superiorization 1181 to implement the regularization reconstruction of 
SPECT. 

The superiorization of iterative methods, which was first proposed by the authors of [16 , is a relaxation 
technology for constrained optimal problem. The superiorized algorithm lies between the feasibility-seeking 
algorithms, which seek a feasible point in the constrained set, and the optimal algorithms, which seek the 
maximum(minimum) point of objective function in the constrained set. The aim of superiorization algorithm 
is to look for a superior point of the objective function instead of the optimal point or just a feasible point 
in the constrained set. The basic idea of superiorization is to do the feasibility-seeking iteration method with 
perturbations from the objective function. 

The superiorization of ART-like algorithms has been studied and applied to the total variation regularization 
reconstruction of computed tomography(CT) |Il]-[2n]- The authors of pLS, Jii, first investigated the convergence 
of the two variants of ART under summable perturbations for consistent case. For inconsistent case, the authors 
of [20j proved the convergence of symmetric version of ART in the presence of summable perturbations of the 
iterates. The superiorization of the EM algorithm was firstly proposed in [21j . and applied to bioluminescence 
tomography. However, to our knowledge, it is still an open problem about the convergence of the superiorized 
EM algorithm. 

In this paper, we will first discuss the convergence of the EM algorithm in the presence of perturbation in 
Section [2j A so-called bounded perturbation resilient (BPR) property of ART is vital in proof of convergence 
for the perturbed version of ART. However, we cannot prove the BPR property of the EM iteration so far, 
because of the nonlinearity of the EM operator. Therefore, we investigate the convergence of the perturbed 
EM algorithm under the following assumptions. Firstly, the perturbations should maintain the positivity of 
iterations. Secondly, the perturbations should go to zero with the increase of iterates. Lastly, the perturbed 
EM iteration should maintain the gradual decrease of the K-L distance between the observed data and the 
projection of estimation by each iteration. 

Based on the convergent conditions, we present the superiorized EM algorithm and its modified version in 
Section [Sj Furthermore, practicable techniques are given to produce ideal perturbations for two special objec- 
tive functions, TV 22 and ^^-norm [53], widely used in imaging science. While the proposed algorithms are 
applicable to diverse inverse problems, in this paper we restrict ourselves to demonstrate its usefulness to the 
SPECT reconstruction problem. Numerical results for SPECT reconstruction are given in Section |4] And as 
expectation, the superiorization algorithms output superior image comparing with the classic EM algorithm in 
terms of MSE and visual quality. Some conclusions and discussions are given in Section [5j 

2 Perturbations Resilience of EM Iteration 

For the sake of reference, wc first introduce some notations and assumptions. In this paper, an image x is 
described as a vector of length N with individual elements xj, j = 1,2, ■ ■ ■ , N. When it is necessary to refer to 
pixels in the context of a 2D image we use the double subscript form Xp^q , where 

J = {q - 1)W + p,p = 1,2,- ■■ ,H,q^ 1,2,- ■■ ,W, (3) 

and integers W and H are, respectively, the width and height of the 2D image array, which has a total number of 
pixels N = W X H . Denote by M.^ the region {x > 0\x S M^}. For the system matrix, we assume Hj — J^i '^ijj 
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and denote by di{x) — ^.ij^j the ith component of Ax for convenience. Furthermore, we introduce 



Let P denote the EM operator, then the EM iteration x''^^ = P{x^) for the problem (|2| is defined as 



with an initial point a;'' > 0. Similarly, the perturbed version of the EM iteration is defined as 



X 



where = x^ + PkV^ ■ Here, the vector v'^ and number /3fc represent the perturbed direction and length for 
the fcth iteration, respectively. Hereafter, we call the EM iteration without perturbation as the classic EM 
iteration. 

It has been established that the sequence {a;"} generated by the classic EM iteration converges to the 
minimizer of the K-L distance l\{x) between b and Ax on i?^, where /^(a;) is defined as 

M , M 

I\{x) = /(6,Ax-)-^6aog— (8) 

where /(•, •) denote the K-L distance function of any two nonnegative vectors. For all vectors a;,?/ > 0, we have 
for /(a;, y) > 0, and the equality hold if and only if a: = y. 

Obviously, if /3fc = 0, the formula ([7| for the perturbed EM iteration is the scheme ([6| for the classic EM 
iteration, in which we are not interested. Therefore, in the following we assume /3fc > 0. A natural question 
is that under what assumptions the sequence a;" generated by the perturbed EM algorithm converges to the 
minimizer of ([8| as well. Before discussing the convergent conditions of the perturbed EM algorithm, we first 
summarize some propositions of the classic EM iteration |H1 [H [23] . Without loss of generality, we assume that 
all the elements hi of h are not zeroes. Furthermore, we assume Hj = Uij = 1 for all j in the theoretical 
derivation for simplicity. 

Theorem 2.1 For any initial point x^ > 0, denote by x'^ the estimate of the EM iteration after k iterations. 
Then the following propositions hold. 

1. x^ > for all k, j, and J^-j = J2i k>Q, 

2. Ia{x^^^) <Ia{x''), and I{x^+^,x^) < I\{x^) I\{x''+^) , 

3. {x^} converges to a minimizer x* of l\{x) on , and x* is a fixed point of the EM operator, 

4. I{x*,x''+^) < I{x\x^). 

The above inequalities hold with equalities in\^ and^ iff is a fixed point of the EM operator. 
Remark 2.2 Obviously, x is a fixed point of the EM operator iff fj{x) ~ 1 for Xj ^ 0. 

From the propositions above, we have that the sequence {/^(a;*^)} monotonically converges to the minimum of 
/^(•) on M^, and {a;'^} approximates to the minimizer x* gradually. We will investigate the convergence of the 
perturbed EM iteration, and prove the similar propositions of the classic EM iteration as far as we can. 

Theorem 2.3 Given any initial point x^ , denote by {a;'^}^,^^ the sequence generated by the perturbed EM 
iteration Q). 



= (4 + (3k ■ v^) ■ f, [x^ + = y'y ■ f, (/) (9) 



where yj—Xj+ (3kVj, and {v'^} are bounded vectors. Several properties of the iterative sequence |^ can be 
identified. 
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1. I\{x^+^) < I\{x^), 

2. {x^} has a convergent subsequence {x™''}, and the limit x is a fixed point of the EM operator, 

3. furthermore, I{x,x^^^) < I{x,x^), and x^ — > x, 

4- lastly, X is the minimizer of l\{x). 

The above inequalities hold with equalities iff x^ is a fixed point of the EM operator. If the following conditions 
hold, 

1. positivity: > 0, 



2. decrease: 
3. 



and 



h max{-^}B- - min{^}B+ + h 



(10) 



where = {j\v^ < 0},S^ = {j\v^ > 0} , and = niax{p + Ejes^ 
min{E^gg+ x^^'^^ , Sjest ^i^' Here p is a sufficiently small positive number. 



Before presenting the proof, we explain the necessity of the conditions of theorem 2.3 The positive condition of 
y'j is necessary for the nonnegative constraints. The second condition is required by the convergence of iterations 
to a feasible point. Intuitively, the last condition is used to guarantee the decrease of the K-L distance l\[x''), 
which means that the iterations should approximate to the minimizer of gradually. 

An important concern is about the exitence of the pertubations satisfying the conditions of theorem |2.3[ 
epecially the third conditio n. B ecause the sequence {x^} generated by the classic EM iteration [fik = 0) satisfies 
the conditions in theorem 2.3 the sequence {x^} generated by the perturbed version will also satisfies them 
when /3fcS are small enough due to the continuities of the K-L dista nce ([s] ) and the EM operator Q. Therefore, 
there exist perturbations which satisfy the assumptions in theorem 1 2. 3 [ 
Proof : We prove this theorem step by step, which follows the proving procedure in [25 . 

1. l\{x^'^^) < I^{x''). Wc will prove that the inequality ( 10 ) implies that {/^(x*^)} is a monontone decreasing 
sequence, i.e. I^{x''^^) < /^(a;*^). By the definition of/^(-), we have 



i 

^ bi log 

i 



d,(yfe)-/?fed.(^;'=) 
di{y 



^+^6,log- 



1 



Mv'')/d^{y'') 



1 - 

diiy'' 



d^jv") 
diivn 



(11) 



> Pk max{- 



_3 



.7GS, 



k 



dijv'') 
d^iy") 



pk max{ 



„.k 



'-k}B-k 



Pk ma^{^^}B^-Pk mm{'^}B+ +PkJ24^'^' 



J6sr 



> Pk max{-^}(B,T - ^ x';+') - Pk min {^}(B+ - ^ 



> /?,max{-^}p-/?, min{^}(B+- ^ x';+') 



jest 



(12) 
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The first and second inequalities hold because of log x > 1— ^ for x > and the condition ( 10 ) , respectively. 
We have proved that J^(x''+^) < l\{x'') by (12 1 and condition 3. Next we should prove that /^(x^^) ~ 



2.1 



l\{x'') iff PkV^ = and x'^ is a fixed point further. Obviously, the sufficiency is true by theorem 
The necessity is also true by the following facts. If l\{x''^^) — l\{x^), we have that all the inequalities in 
derivation above hold with equalities, which imply that i_^^rf.(Jfc)/(j.(t,fc) — 1 and 5^7 = by inequalities 
Vn and (12), respectively, i.e. I5kdi{v^) = and > 0. In fact, we have = 0(V j), i.e. y'^ = x'^, and 



that x'' is a fixed point by theorem 2.1 If there exists jo such that wj^^ > 0, there must exist ig such that 

aigjg > and di^{v^) ~ o.ioj'v^ ^ '^iojo'^jo > by the assumption = 1 7^ 0, which is contradict 

to d,{v^) = 0(V i). 

2. Existing convergent subsequence {a;™'=}: Because < — X^i ('^O'^^tant) and Xj > 0, {x''} 
has a convergent subsequence {x'"'-'}. Denoting by x the limit of {a;™'=}, we have y™*" — > x because 
/3fc — > and v'^ is bounded. 

Next we will prove i is a fixed point of EM operator. To this end, we define a function about a; > 

D{x) = I{P{x),x). (13) 

By the second proposition of theorem |2.1[ we have 

Diy"^") = /(P(y"'=),y'"'-) - /(x™''+i, < - /^Ix^^+i). (14) 

Due to the statement 1, {/^(x'')} is a convergent sequence. Furthermore, l\{y^) converges to the 
same limit as l\{x^) because of the continuity of and (5k — > 0. Therefore, we have l\{y'^'') — 

l\{x"^>'+'^) — > 0, TOfc — ^ oo. Thus, we have that 

D{x) = lim Diy""-) = I{P{x),x) = 0, (15) 

nifc >-oo 

i.e. i is a fixed point of the EM operator, since I{x,y) = ■^==> x = y. Therefore, we have fj{x) = 1 if 
Xj 7^ by the EM iteration formula. 

3. /(x,x'^+i) < I{x,x''). Before going further, we first prove a general inequality I{x,x) — I{x,P{x)) > 
l\{x) — l\{x) for any point x > 0. Assume z = P{x), then we have 

/(x, x) — /(x, z) 

Xj\0g^~ V(Xj - Xj) 



3 3 



/ J "^j / ^ y^3 / \ / - \ / . 



Z3 9tjix) 9ij{x 



5,,(x) log - V(x, - X,) 



^ X, ^ g., (X) log g 1^ + E E ^'^^ ) if) - E 



> - E X, log E (x) ^ 1^ + E E (^) log " E - ) 

X, log ^L^^^ + \ \ Xj -y— log — ^ - ) [Xj - Xj) 



3 



» 3 

i ' i ' i i j 

= I\{x)-I\{x), (16) 

where we used % ~ ~ Tli^i ~ ^'^'^ the definition of /(•,•). The inequality and the 

second equality were implied by Jensen's inequality and the fact fj{x) — '^^gijix) ~ 1 for Xj 7^ 0, 
respectively. The last equality follows J2j^3 — ^^'^ J2idiix) — '^ij^j ~ the 



assumption J2i^ij — ^ ^^'^ the property 1 of theorem 2.1 Let x — y'^, then z — P{y'') — x''+^ and (16) 
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implies 



I{x,y'^)-I{x,x''+') > I\{y^)-I\{x) 

= [I\{V^) - l\{x'^+')\ + I\{x^+') I\{x) 

> I\{y^)-I\{x'^+^). (17) 

The last inequality held because l\{x^'^^) > I^ix)- 

Next we have prove the property of {x*^} that I{x, x''^^) < I{x, x^). By the definitions of the K-L distance 
and the perturbed EM formula, we have 

I{x,x^)^ I{x,x^+^) 

{x^ + Pkv^)fj{y'') 



x^ 



- J2 S:, log /, (/■) + J2 ^1 + /^fc ^ ) 

X A 

3 J ^ 

> yi iogyHi^ + YiJi J^^] (18) 

= /(.i,,'=)-/(i,.^+i) + ^f,.^^_/3,^.| 

> 4(/) - /^(x'^+i) - maxj-"^} V X, + /3, miii {"^j V x, - V z;}' 



> niax{-4:}(i?- - ^ i,) + fS^ min {4:}( ^ - B+) 

> 0, (19) 



Again we used the inequality logx > 1 — - for a; > for the first inequality. The second and the third 
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inequalities hold by (17) and amplification and shrink techniuqe. The last inequatility hold by the formula 



We have that I{x,x^'^^) < I{x,y^) by equation (19 1. Furthermore, we have /(i, 2;''+^) = I{x,y'') iff 



Vj=0 and x'' is the fixed point of the EM operator. The sufficiency is clear. 
The necessity is also true by the following facts. If I{x, x'^+^) = I{x, y*^), we have that all the inequalities 
in deviation above hold with equalities, which imply that w| = and l\{y^) = l\{x^'^^) by inequalities 



2.1 



( 18 1 and ( 19 ), respectively. Therefore y = x and x is the fixed point of the EM operator by the theorem 



Now we prove the convergence of {x^} generated by the perturbed EM formula, i.e x^ — > x. F rom 



(19 1, we can conclude that I{x,x ) is a monotone decreasing sequence if the assumptions of theorem 2.3 
hold. Thus, we have I{x, x^) \ because of /(i, x™^ ) — > and the monotonicity of /(i:, a;''). Therefore, 
we have proved x^ — > x, and y'^ — > x by /3fcw'^ — > and I{x, y) — x = y. 

X is a minimizer of on M" . In order to prove this statement, we should prove x satisfies the Kuhn- 
Tucker (K-T) conditions. For the EM algorithm and the K-L function, the K-T conditions are equivalent 
to (1) fj(x) = 1 ff ^ and (2) < fj{x) < 1 if = P Hi]. Obviously, the firstly condition holds 
because x is the fixed point of the EM operator. For the second condition, the nonnegativity of x is 
satisfied since x is the limit of x'' > 0. Next we prove that fj{x) < 1 for Xj — 0. If fj{x) > 1 + e > 1 
for Xj = 0, by the continuities of fj{-) and the EM operator and y'' — > x, wc have that there exists a 
sufficiently large integer L > such that fj{y*') > 1 + ei with ei > for all A: > L. If x^ — > 0, there must 
be infinite iterations such that x^"*"^ ^ ^j- ^^'^ following, we assume that k > L. For each iteration 
satisfying x'^^^ — y^ ■ fj{y^) < Xj, we have that 



1 > 



= (l + /3/c4)/.(2/') (20) 
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By the assumption fj{y^) > (1 + ei), we have 



l+/3fc^<T— ■ (21) 



Abstracting 1 on both sides of (21), we have Pk:^ < — jq:^- By equation (12 1, we have that 



> (22) 
1 + ei 

The second inequahty hold because yj=Xj+ < < and < J2jes+ ^i^^- "^^^ ^^^^ inequality 
follows from the definition of Bj^ . 

Tb l^k\ 



From equation (22), we can see that the decrease of I^{x ) is larger than a positive number for each 



iteration satisfying x- < Xj. Further, because the number of iterations satisfying x^ < Xj is infinite 
and l\{x^) is finite number, /^(a;') will go to negative infinity from, which is contradict to the fact that 
I\{x^) > for all xK 

3 Superiorization of the EM Algorithm 

The regularization reconstruction for SPECT is often modeled as a constrained optimization problem 

min(^(a;), £■ = {x*|a;* argmin/^(a;)}, (23) 

xGE x>0 

where <^ is a convex function, which assigns an image x a number indicating the " undesirability" of the image 
in some sense. The set E is called feasible set, and it is called feasible problem to find a point in E. 

To our knowledge, there is no efhcient algorithm to deal with the constrained optimization problem ([23]) 
because of the large scale of it for SPECT reconstruction. On the other hand, although the feasible problem 
is also a constrained optimization, we can solve it by the classic EM algorithm [5] and its variants (TUl E] 
efficiently. Based on the facts above, we use the superiorization methodology for the given objective function 
to implement the regularization reconstruction of SPECT. 

For the objective function 0, the superiorized EM algorithm is illustrated as algorithm [T] based on the 
conditions of theorem |2.3[ In order to emphasize the objective function for which we are superiorizing, we 
refer to the superiorized algorithm as the 0-superiorization of the EM iteration. 

Algorithm 1 Framework of 0-supcriorization algorithm 

Initialization: /3o > 0, a;° > 0, fc = 0, and < 7 < 1. 
repeat: logic=true 
while logic 

find a decreasing direction u*^ of cj) at x*"', such that yj—Vj+ jikv^ > 0; 



If (jjiy'^) < (/'(a;'') and inequality (10) hold (*) 

logic=false, x''+^ = P{y'"), Pk+i ^ Pk, k ■ 
end(if) 

Pk = llik 
end(while) 



Because the perturbation direction v is selected as the decreasing direction of (p at x , the size of /3fc 
represents the strength of regularization in some sense. Because the condition 3 of theorem |2.3| is very strict, 
the numerical experiments show that {/3/c} goes to zero very fast, which results the the regularization effect is 
_very weak. Therefore, we propose a modified version of algorithm[l which is show n in algorithm [2| In algori thm 



[2] we only validate 



k + l 



implies l\{x^'^'^) < l\{x^) by theorem 



2.3 



i[Tl wl 



) < /^(a;'^), rather than the inequality {Um of theorem 2.3 Since the inequality (10) 



algorithm 2 can be seen as a relaxation of algorithm [l] Furthermore 
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we introduce a relative decrease of to avoid the situation that the amount of /^(a;*^) — /^(x'^^^) for each 
iteration is too smaU, which can accelerate the convergence of l\{x^) and x'^ intuitively. 



In order to confirm the inequality (10) of theorem 2.3 in algorithm 111 we should compute X^jesT ■^i ^"^^ 



Algorithm 2 Modification of the 0-superiorization algorithm [T] 



/3o > 0, a;° > 0, fc = 0, and < 7 < 1. 
repeat: logic=true 
while logic 

find a decreasing direction of (p at x^ , such that 2/1="^^+ PkV^ > 0; 
If < (l){x^) and I'XiPy^) < Ia{x^) W 

logic=false, = P{y^)- 

Pk+i = iPk (02) 

else 

— Pk 

end(else) 
fc = fc + 1 
else 

Pk = 7/3fc (03) 
end(else) 
end (while) 

Sjgs+ fo'^ each iteration. However, the limit x is not known in the computational procedure. In practice, 
we estimate the values J^j^s^ ^^"^ Sjes+ % by 



E - (24) 



E - (25) 
where jS*^! and jS"^! denote the cardinalities of the sets and , respectively. Furthermore, the capital 



letters TV and B = formulae (24 1 and (25) denote the length of x and the total counts, respectively. 

Further, the parameter p in theorem 2.3 is chosen as 0.01 and 10~^. 

In the following, we will discuss how to generate desirable perturbation f3kv'' or y'^ for two concrete ob- 
jective functions, TV and /^-norm, such that the perturbations satisfy the conditions (*) and (*) of the two 
0-superiorization algorithms. Since the TV regularization allows the reconstructed image to have sharp edges, 
TV based models are widely used in imaging sciences [THl [221 123 HZ] ■ For an x T4^ image x whose pixel values 
are denoted by Xij, the TV of x is defined as 

H-lW-l 



TV{x) = ^ ^ Y'(xi+ij-a;,j)2 + (xij+i-a;,j)2. (26) 

where H, W are the height and width of x. In order to decrease the value of TV at x'', we choose v'^ as 

i;'= = sV|s1 (27) 

where s'' £ dTV{x*') is the sub-gradient of TV at point x^ , and \s^\ is the maximum absolute value of the 
components of . Therefore, the sequence {v^} is bounded. In fact, is the normalization of in the 
space, rather than in the P space used in [T8l [20] . 

In addition to the TV based models, Z^-norm minimization method is another widely used technique in image 
sciences and processing [IS] [25] US] • Here the /^-norm is about the wavelet coefficients of a;, which is defined as 

\\T{^M^)\\i = El"^l (28) 
j 

where UjS are the coefficients of x under a given wavelet basis ^nd the letter T denotes the wavelet 

decomposition operator. Although we can use the same method for TV function to reduce the wavelet ^^-norm 
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(a) (b) 
Figure 1: (a): activity map (b): attenuation map. 



of x'^ , we introduce two more effective methods, soft and hard thresholding schemes, to decrease the wavelet 
l^-norm of x'^ . Define 7^ by 

f \a^\ > pk 

7t = Hard(a,^) = ' , (29) 



. -sign(a*^) la^l > /S^ 
7? ^ S„ft,a?,H 51 ,4,^,^ ■ (30) 



where a'j are the wavelet coefficients of x''. The perturbation direction for each iteration is selected as v'^ = 
T^^y{'^^), and y*^ = a;'"' + /Skv'^. In fact, we need not to compute v'^ explicitly. By the linearity of wavelet 
transform, we can compute the wavelet coefficients of y'" directly by 

a; - H„d-(af).{°J , (31) 

aj = S„f.-(a5) = |°J-='P'<»5-)-A , ,32) 

and — T^^ y{a'j). In the numerical simulations, we use Daubechies 6.8 biorthogonal wavelets with symmetric 
extensions at the boundaries. For refereed convenience, we use hard-superiorization and soft-superiorization to 
distinguish them for Z^-superiorization in this paper. In order to avoid < 0, y^ is set as ^x'j if < in 
practice. 



4 Numerical Results 

In this section we investigate the stability and robust of the proposed algorithms by several numerical experi- 
ments of SPECT reconstruction. To this end, we conducted simulation studies by employing computer, and the 
projection data were generated based on the following model. Figure [T] (a) and (b) display a chest phantom of 
activity and an attenuation map, respectively. The chest phantom activity consists of an ellipsoidal background 
(body region) with axes of length 22.5cm and 30cm, which contains two smaller ellipsoidal regions (lungs) with 
axes of length 10cm and 8.8cm, and a ring(myocardium) which inner and outer diameters are 6cm and 8cm, 
respectively. The activity in myocardium, background, and lungs were specified to be in the ratio 8:3:1. 

To simulate the strongly non- uniform attenuation coefficient of chest we utilize the phantom used in [TU], 
which imitates the attenuation map across a section of human thorax. Besides the body region and the lungs 
of the activity map, the attenuation map consists of two circular regions (bones) of diameter 2.5cm(see figure 
1(b)). Attenuation coefficients were 0.03cm~^ within 'lung' regions, 0.17cm^^ within 'bone' regions, 0.15cm~-'^ 
elsewhere within the body ellipse, and O.OOcm"^ outside the body. 

In the simulation, the activity and attenuation maps are evenly sampled in [—15, i5] x [—15, 15] on a grid 
of 128 X 128. A perfect parallel hole collimator was assumed, and noise-free projection data were created via 
attenuated Radon transform formula ([I]), which included tissue attenuation, but scattering and blurring. In 
order to simulate the quantum noise in the simulated data, the following procedure was implemented for each 
data set [3T] . 
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Figure 2: Images reconstructed from the simulated data set l(top row) and 2 (bottom) by the classic EM 
algorithm(left column) and the mean images of 100 trials of noise(right column). 



• The projections were scaled (multiplied by a constant factor) so that the total count was a predefined 
integer, 

• each point value in the data set was then replaced by a random realization of a Poisson variate with a 
mean equal to that value. 

Two data sets, sixty and thirty projections, were generated over 180° evenly with view angles (pi = ^^7r(Z — 
1, • • • , A^'o, and A^o ~ 30 or 60). Counts were recorded in 128 bins per projection. Total counts recorded on all 
projections are about 500K and lOOK for two simulation data sets. 

In order to evaluated the performance of the superiorized algorithms, we computed the mean image x* 
from 100 noise trials of the two data sets. 

^ 100 

X* = y i", (33) 

100 ^ ^ ^ 

m— 1 

where x"^{i = 1, 2, • • • , 100) is the reconstructed image of the mth noise trial by the classic EM algorithm, and 
the number of iteration for each trial is 30. 

Figure [2] displays the reconstructed images by the classic EM algorithm and the mean images for two 
simulated data sets. As expected, the mean images for the two data sets are clear visually comparing with 
the images reconstructed by the classic EM algorithm, which are dominated by noise, especially the image for 
data set 2 due to the lower counts. Thus, we used the mean square error(MSE) between the estimations and 
the mean images for two data sets to measure the image qualities. And we used it as the criteria to to select 
the outputs of the iterative algorithms. The mean square error(MSE) between the estimation x and the mean 
estimation is computed by 



1 ^ 

MSlL{x) ^ -Y.^x, ~ x*f (34) 



The outputs of the different algorithms are taken as the bets estimations in terms of MSE. In order to compare 
the image quality from different data set, we introduce an amount of relative MSE(RMSE), which is defined as 



RMSE{x) = 



where x* is the mean image of the same data set. 

In the numerical simulations, an initial image with value c ~ hi is used for all experiments, unless 
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there is a further explanation. The parameters /3o are chosen as c/2 and c/10 for the TV- and Z-'^-superiorized 
EM algorithms, respectively. The threshold Hi of the relative decrease of /^(a;*^) in the superiorized algorithm 
[2] and the parameter 7 are chosen as 0.01 and 1/2, respectively. 
Experiment 1: sixty projections and 500K counts 

Figure [3] displays the images reconstructed from data set 1 by the TV-, hard- and soft-superiorized EM algo- 
rithms from column 1 to column 3, respectively. Figured plots the variances of MSE(a:;''), /^(a;*^) and Pk versus 
the iteration number. Table [l] tabulates the TV values, 7 -norms, RMSEs and iterations of the images in figure 
[3j the image by the classic EM algorithm and the mean image for data set 1 in figure [2] In the following, 
TV-, hard- and soft-n(n — 1,2) are the abbreviations of the TV-, hard and soft-superiorized EM algorithm 
n(n = 1, 2) for simplicity. 





Figure 3: Images reconstructed from data set 1. The images in top and bottom rows are reconstructed by 
the superiorized EM algorithm [T] and algorithm [2] respectively. The images reconstructed by TV-, hard- and 
soft-superiorized EM algorithms are displayed from column 1 to column 3. 



Table 1: TV, Z^-norm, RMSE and iteration of the images in figure 



image 


EM 


TV-algl 


hard-algl 


soft-algl 


TV(xlO^) 

RMSE 
iterations 


26.708 
12.953 
0.1914 
13 


25.318 
12.404 
0.1873 
13 


25.577 
12.320 
0.1875 
13 


25.958 
12.489 
0.1888 
13 


image 


mean 


TV-alg2 


hard-alg2 


soft-alg2 


TV(xlO^) 

RMSE 
iterations 


13.987 
8.5244 


19.434 
10.279 
0.1633 
15 


14.115 
2.867 
0.1563 
30 


15.443 
3.189 
0.1747 
23 



As expectation, we can see that the images by the superiorized EM algorithms are superior to the one 
by the classic EM algorithm from the comparison of TV, Z^-norm and RMSE values in table [ij though the 
images by the superiorized algorithm 1 are visually indistinguishable from the one by the classic EM algorithm. 
Furthermore, the ^-superiorized EM algorithm[2]are superior to the corresponding ^-superiorized EM algorithm 
1. And we can see that the results by the Z^-superiorized EM algorithms [2] are superior to the TV-superiorized 
algorithm [2] visually. The reason for above phenomena maybe that the parameter /S^ goes to zero very fast for 
the superiorized EM algorithm [l] and the TV-superiorized EM algorithms 2 (see figure |4]) 



From figure |4| we can see that the l\{x'') decreased if Pk and w'"' satisfied the conditions 3 of theorem 



which validate the conclusions of the theorem 



2.3 



2.3 



numerically. Further, we can see the plots of /^(a; ) for 



the classic EM algorithm and the superiorized EM algorithms [T] are indistinguishable from each other because 
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Figure 4: Plots of the K-L distance /^(a;'^)(top row), mean square error MSE(2:'') (middle row) and /3fc (bottom 
row) versus the iterations k of the classic EM iteration and the superiorized versions(left for algorithm 1 and 
right for algorithm 2). 



the parameters /3fe went to zero very fast for the superiorized algorithm 1. The reason for the fast decreasing 
of l3k for the superiorized algorithm 1 is that the inequality for condition 3 is very strict because of the use 
of magnification technique several times. Lastly, we can see that the plots of MSE(a:;'^) for the classic EM 
algorithm and the superiorized version [l] are presented as L-curve: decrease first and then increase, though 
l\{x'') decreased gradually. And we can see that the plots of MSE(a;'^) for the -superiorized EM algorithm 2 
decreased unceasingly, which implies the estimate x'^ for each iteration approximate the mean image gradually. 

From table [l] we can see that the iterations of the superiorized EM algorithm 2 is lager than the classic 
EM and the superiorized EM algorithm 1. However, the RMSE of the estimation of the TV-, hard- and soft- 
superiorized algorithms 2 after 13 iterations are 0.1668 0.1741 0.1797, respectively. Therefore, the superiorized 
algorithms 2 are superior to the classic EM and the superiorized EM algorithms. 
Experiment 2: thirty projections and lOOK counts. 

Figure [5] displays the images reconstructed from data set 2 by the TV-, hard- and soft-superiorized EM 
algorithms 1 and 2 from column 1 to column 3 and from row 1 to row 2. Figure [6] plots the variances of 
MSE(a::'^), l\{x'') and Pk versus the iteration number. Table [2] tabulates the TV values, ^^-norms, RMSEs and 
iterations of the images in figure [5l the image by the classic EM algorithm and the mean image for data set 
2 (see figure [2]). 
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Figure 5: Images reconstructed from data set 2. The images reconstructed by the TV-, hard- and soft- 
superiorized EM algorithm 1 and 2 are displayed in column 2, 3 and 4, and row 1 and 2, respectively. 



Table 2: TV, l^-norm, RMSE and iteration of the images in figure [sj 



image 


EM 


TV-algl 


hard-algl 


soft-algl 


TV(xlO^) 


13.739 


13.105 


13.396 


13.061 


6.729 


6.452 


6.451 


6.277 


RMSE 


0.2822 


0.2778 


0.2802 


0.2795 


iterations 


8 


8 


8 


8 


image 


mean 


TV-alg2 


hard-alg2 


soft-alg2 


TV(xl03) 


7.167 


9.769 


5.733 


5.576 


4.037 


5.069 


1.434 


1.283 


RMSE 





0.2490 


0.2091 


0.2064 


iterations 




9 


12 


10 



From the observations of figure [5]and|6l and table [2) we can draw the same conclusions as the experiment 
1. Further, we can see that the visual quality are inferior to the images of experiment 1 by comparing with 
figure [3] because of the low counts level. And the superiority of the (^-superiorized EM algorithm 2, especially 
the ^^-superiorized algorithm 2, is much more remarkable than the experiment 1. 

For the two experiments above, though the Z^-superiorized EM algorithms obtain better results than the 
classic EM algorithm and the TV-superiorized versions, the reconstructed images by Z^-superiorized EM algo- 
rithms contains Gibbs phenomenon caused by the hard- and soft-threshold techniques. 

From figure [2] [3] and [5j we can see that the right part of the images are blurred strongly. This is because the 
detectors rotated from up to bottom through right, and the the gamma rays emitted from the right hand pixels 
were more likely to be absorbed in propagation. 

In the following experiment 3, we will see the outstanding performance of the superiorized algorithm 2 for 
initial vector with random elements. 

Experiment 3: initial 2:° with random values on interval [1, 2] for data set 1 

This experiment is used to test the stability and robustness of the superiorized EM algorithms for initial 
point x'^ . In this experiment, the initial point x'^ is selected as a random vector with values on interval [1,2]. 
The reconstructed images by different algorithms are displayed in figure [Tj and the corresponding TV values, 
Z^-norms and RMSEs are tabulated in table H] 

We observe that the TV- and Z^-superiorizcd EM algorithms 2 are stable and robust for initial point. By 
comparing the figure [3] andjTj and tabic [T] and |3] we have that the outputs of the classic EM algorithm and 
the superiorized algorithm [ijare strongly effected by the initial points, while the effect of initial point for the 
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Figure 6: The K-L distance l\{x''){top row), mean square error MSB (x*^) (middle row) and /3fe(bottom row) 
variance versus the iterations k of the classic EM iteration and the superiorized versions (left for algorithm 1 
and right for algorithm 2) for experiment. 



TV- and Z^-superiorized EM algorithm 2 is very weak. Because the variances of l\{x'') and MSE(x'') versus 
iteration are very similar to that of experiment 1, we only plot the variance of /3fc in figure [8] 

A surprising observation is that the random initial vector is superior to the uniformly initial vector for 
the superiorized EM algorithms 2 in term of RMSE by comparing the table [l] and |3] This is contrary to the 
long-standing initial selection of iterative method for image reconstruction. 

By comparing the plots of /3fc of the experiments above for different algorithms, we can observe the following 
facts. Firstly, the parameter /S^ went to zero very fast for the superiorized EM algorithms 1. Secondly, the 
parameters of experiment 1 and 2 went to zero slowly for the /^-superiorized algorithm 2 while fast for the 
TV-superiorized EM algorithm 2 with uniformly initial point (see figure |4] and |6]) . Thirdly, the parameter 
of experiment 3 went to zero at a relatively slow rate for the TV-superiorized EM algorithm 2 comparing with 
experiment 1 and 2. 

The reasons for this include the following aspects. Firstly, as discussed above, the parameter /3k for the 
superiorized EM algorithm 1 went to zero fast by the strict condition (3) of theorem 2.3 Secondly, the estimation 
of each iteration with uniformly initial vector x° is very smooth, which results the parameter f3k goes to zero 
fast by the decrease condition of the TV function in formulae (*) and {-k). On the other hand, the estimation by 
the TV-superiorized algorithm is very rough, which results the decrease of the TV function for relatively large 
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Table 3: TV, wavelet ^^-norm, RMSE and iteration of the images in figure [sj 



image 


EM 


TV-algl 


hard-algl 


soft-algl 


TV(xlO^) 

RMSE 
iterations 


40.591 
19.098 
0.2538 
13 


29.804 
14.632 
0.2112 
13 


40.364 
18.907 
0.2529 
13 


34.214 
16.138 
0.2266 
13 


image 




i V-alg2 


1 1 „ 1 „o 

nara-alg2 


sott-algi 


TV(xl03) 

RMSE 
iterations 




17.133 
9.403 
0.1469 
19 


16.203 
3.417 
0.1683 
30 


16.372 
3.207 
0.1873 
23 



parameter /3k- Thirdly, the threshold techniques (hard and soft) always decrease the ^^-norm, which results the 
conditions (*) and (★) for the ^^-superiorized algorithm 2 become to only one condition about the decrease of 
K-L distance. 

Experiment 4: modification of the TV-superiorized algorithm 1 and 2 

This experiment is based on the data set 1. We modify the the TV-superiorized algorithms by discarding the 
deceasing condition of TV in formulae (*) and (*). 

Figure [9] displays the reconstructed images by TV-superiorized EM algorithm 1 and 2 in absence of the 
decreasing condition of TV function. The quantitative values are tabulated in table |4j And figure 10 plots the 
variances of Pk versus iteration number. 

As expectation, the parameter f3k converged to zero at much slower rate comparing with the TV-superiorized 
algorithm 2, which is similar to the ^^-superiorized EM algorithm 2. However, this modification have very little 
effect on the TV-superiorized algorithm 1. The observations above show that the decreasing condition of TV is 
dominated for the TV-superiorized EM algorithm 2, while the condition 3 of theorem 2.3 for the TV-superiorized 
EM algorithm 1. 

It is amazing that the modified algorithm also decrease the TV function, even the reconstructed image is 
better than that the TV-superiorized EM algorithm 2 , though we do not validate the decreasing condition of 
it. 



Table 4: TV, wavelet Z -norm, RMSE and iteration of the images in figure §(xl03) 



image 


image 1 


image 2 


image 3 


imagc4 


TV(xlO^) 


25.736 


14.101 


13.105 


7.758 


/i(xl03) 


12.534 


7.494 


6.4517 


3.6920 


RMSE 


0.1884 


0.1278 


0.2778 


0.2075 


iteration 


13 


27 


8 


16 



Remark 4.1 The experiments show that the convergence of the superiorized EM algorithm 2, though we can 
not prove it theoretically so far. 

Remark 4.2 The parameter jS^ represents the strength regularization of (j), so we can obtain the regularization 
reconstruction by terminating the iteration as long as {3^ is smaller than a predefined number e. The discussion 
about the selection of e will be studied in future work. 

5 Conclusions and Discussions 

In this paper, we discussed the convergence of the EM algorithm in the presence of perturbations, and proposed 
the superiorized EM algorithm [l] based on the convergent conditions. The modified version of algorithm [l] was 
proposed by relaxing the condition 3 of theorem |2.3[ We gave the methods to generate the desired perturbation 
for TV- and /""^-superiorization algorithms. The numerical experiments for SPECT reconstruction show that the 
superiorization algorithms is more stable and robust than the classic EM algorithm for low-count data set and 
initial point, and they could efficiently decrease the corresponding objective functions as expectation. 

A lot work should be developed to about the superiorized EM iteration. Could we prove the convergence of 
the superiorized EM algorithm 2 theoretically? Secondly, we have not exact theories to support the amazing 
observation of the experiment 4. The observations of the numerical experiments, especially experiment 3 and 
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Figure 7: Images for experiment 3. The top image is reconstructed by the classic EM algorithm. The images in 
the column 1, 2 and 3 and the middle and bottom rows are reconstructed by the TV-, hard- and soft-superiorized 
EM algorithm 1 and 2, respectively. 

4, would enlighten us to modify the the superiorized of the EM iteration. 

In addition, the concept of superiorization is an emerging technique, and a lot of works need to be done 
and developed. For instance, we could not prove (t>{x*) < 4>{x) theoretically, where x and x* are the solutions 
by the classic iteration algorithm and the (/)-superiorization version |18) . Secondly, the superiorization theories 
are mainly about algebraic iteration [T5| so far, and it is difficult to generalize the theory about algebraic 
iteration to other iteration algorithm. Specially, due to the nonlinearity of EM algorithm, we can not prove the 
bounded perturbation resilience(BPR) of EM algorithm, and the convergent conditions in this paper for the 
perturbed EM algorithm are much stronger than that for the perturbed algebra iteration. 
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